THEORETICAL APPROACH TO THE CORRELATIONS OF A HARD-SPHERE 
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We present the first theoretical approach to the angular- 
average of the two-body correlation function g(r) for simple 
solids. It is based on three sum rules for g(r): the compress- 
ibility and virial equations and the normalization. We apply 
the theory to determine this correlation function for the case 
of the FCC solid phase of hard spheres. The agreement with 
simulation data is excellent over all the density range. The ap- 
plication to other simple systems is discussed. The approach 
opens a new route to perturbation theories for simple solids. 

PACS numbers: 61.50.Ah,64.70.Kb 



The most important correlation function in fluid clas- 
sical systems is the radial distribution function g(r). The 
theoretical accessibility of this function for a hard-sphere 
(HS) fluid, namely the Percus-Yevick approximation, has 
been the crucial point to the extraordinary development 
of the theory of simple liquids 0^). The analogous 
progress is lacking in the case of structured phases of 
these systems, the solids. Spite of the advance in the the- 
ory of classical non-uniform systems experimented dur- 
ing the last fifteen years ||, one of the most important 
theoretical objectives, correlations, has remained unap- 
proachable. 

In structured systems the most important correlation 
is g(r) the angular average of the two-particle density 
function p^ 2 \r, r'). It appears in all crucial equilibrium 
equations as the virial and compressibility equations or 



the energy equation 
tion approaches J 



and it is the key of all perturba- 
The g(r) plays the same rule in 
classical solids as g(r) does in classical fluids. However, 
up to now, no theoretical approach to g(r) has been re- 
ported. Here, we present the first theoretical approach 
to this function and apply it to evaluate g(r) for a face- 
centred-cubic HS solid. The agreement with simulation 
results is excellent over all the physical range of densities 
from below melting up to almost close packing. 

All previous theoretical approaches reported in the 
literature to the correlations of classical solids have 
been focused on g(r,r') (defined through p^ 2 \v,v') = 
p(r)p(r')g(r,r')) for HS and mainly in relation to per- 
turbation schemes. These are not true approximations, 
rather, they are simple mappings to the radial distribu- 
tion function of a uniform system at some effective den- 
sity. None of these approaches can resist a direct compar- 
ison with simulation results H. In fact, these mappings 



could be done over a wide family of functions || . More- 
over, only p( 2 )(r, r') (the probability to find two particles 
at r and r') and p( 2 \r,r')/ p(r) (the probability to find 
a particle at r provided there is another one at r') and 
their angular averages have a direct physical meaning. 
Then it would be more sensible to make approximations 
on these functions than on g(r,r'). The function g(r) is 
defined as 
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where V is the volume, p = N/V is the mean density and 
dCl the differential solid angle aperture around r. 

To motivate the new theory we shall first discuss the 
two types of correlations that g(r) must accounts for in 
a solid. The long-range correlation characteristic of the 
structured phases is directly due to the periodicity of the 
solid. If no other correlation is considered each particle 
would move around a lattice site independently of each 
others. The probability of finding a particle at r is given 
by the local density p(r) which is a sum of Gaussian-like 
functions located at each lattice site. Then p( 2 \r,r') is 
given by the simple product of individual probabilities 
p(y)p(y') multiplied by a step-like function, g(r, r'), to 
avoid the double occupancy. Thus, the long-range corre- 
lation is already described by the product p(r)p(r'). Its 
proper angular average gives the long-range contribution 
to g(r). Then, it is convenient to define <7o( r ) as the 
angular average of the product p(r)p(r'), 



9o (r) = 
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dfl / dr 1 p(r 1 )p(r 1 + r) 



(2) 



which has the form of a sum of Gaussian-like peaks g^ 
centred around successive neighbours distances Ri, nor- 
malized to the corresponding number of neighbours nt. 
Then, if no other correlation than the long-range correla- 
tions imposed by the lattice is considered, g^> — g^ for 
i > 0, and both functions are the same except that the 
first peak of go(r) is excluded in g(r). This long- range 
correlation is the only one existing in the close-packing 
limit of a HS solid. In this limit p^ 2 \v,v') is exactly a 
sum of products of 5 functions at any pair of different 
lattice sites multiply by a step-like function to exclude 
the double occupancy. Hence, the function g(r) will be 
a sum of 8 functions at different neighbours distances, 
except at zero distance, normalized to the corresponding 
number of neighbours. 



To take into account other correlations the successive 
peaks must be increasingly modified from the g^ 
i decreases to 1. The compressibility equation [0 



as 



Anp / drr [g(r) - g (r)} = -1 + pk B T\T- 



(3) 



is a sum rule for g(r) which help us to understand these 
modifications. The integral of go(r) in Eqn.(^|) gives the 
number of particles inside the system provided the origin 
is located at a lattice site. The integral of g(r) gives the 
number of particles inside the system (minus one because 
of the self-exclusion) provided a particle is fixed at the 
origin. Thus the left hand side of Eqn. (|J) can be under- 
stood as the number AN of particles (minus one) com- 
ing into the system when a particle is fixed at the origin, 
therefore AN = pk B T\T. Observe that in the HS close- 
packing limit xt = 0, the correlation reduces to the long- 
range one discussed above, hence = g^' and Eqn.(|^) 
is verified identically. In other words, no HS come into 
the system if it is completely packed. If we imagine a 
spherical system of radius R we can estimate AN as 
p4:irR 2 5R, where 8R is the displacement of the peaks at 
the border of the system, i.e. the displacement of g^ re- 
spect to g^ for R t w R. Then SR « pk B T X T I {p^Rj ). 
This is a quite interesting result which shows that the 
differences, which should include some kind of deforma- 
tion, between and g^ reduce quadratically with the 
distance. 

Until here the discussion has been quite general and 
now we apply it to the case of a HS solid. The com- 
pressibility of the HS solid is so small that 5R would be 
practically imperceptible in a simulation. A rough esti- 
mation of this displacement is already negligible for the 
first peak even at the lowest densities: taking R w 1 and 
pk B T X T ~ 0.02 give 5R « 0.001 (distances in HS diam- 
eter units dns)- Thus the location of the peaks of g(r) 
does not differ from the localization of those of go (r) . If 
the displacement of the peaks is negligible it is quite sen- 
sible to assume that their functional form cannot differ 
significatly from the peaks of go{r) except for the first 
one due to the characteristic cutoff of the HS. 

Up to here we have described the subjacent physics 
which motivates the following theoretical approach to 
g(r). Accordingly we approximate the peaks of g{r) by 
those of go(r) beyond nearest- neighbours. To go far to- 
ward an explicit form for these peaks, we can regard the 
existent functional theories |s| which give an accurate de- 
scription of the free energy of a HS solid and from which 
it is possible to determine p(r). All these theories use 
successfully the Gaussian parametrization of p(r) 0, 
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where a is the Gaussian width parameter. Using Eqn. 
into Eqn.(0) yields: 



-a(r-Ri) 2 /2 



-a(r+R.i) 2 /2 



rR, 



i > 0. 



(5) 



For the sake of simplicity, we have dropped all terms 
which arise from the exponential products with R R' 
in (||). At the usual values of a they give negligible con- 
tributions because of the absence of overlapping. 

The above discussion on the compressibility suggests 
for g^(r) a functional form similar to that of Eqn.(|^). 
We propose the simple parametric form: 



~gW(r) = 



A e -ai(r- ri ) 2 /2 



> d 



HS, 



(6) 



with g^'(r) = for r < d B s- The compressibility 
shows that even in the less favourable case, at the low- 
est densities, the displacement of the first peak is quite 
small. This suggests that mean location of the nearest- 
neighbours < r > can be approximated by the mean 
value obtained with the first peak of go{r): 



< r >= I drrg^ir) = 



dvr~g { 1] {r). (7) 



This equation is a sum rule for g(r) becomes less approx- 
imate as the mean density increases and it is exact at the 
limit of close-packing. Besides two other exact sum rules 
must be obeyed by g(r). The first one corresponds to the 
normalization of g^ to the nearest-neighbours number: 



47T/9 



drr 2 g^\r) = n±. 



(8) 



dns 



The virial equation is the second exact sum rule. It can 
be easily proved that, for non-uniform systems, the pres- 
sure is related to the value of g(r) at contact exactly in 
the same way as it is related to the radial distribution 
function of uniform fluids: 



f3P/p=l + 47 1 g(d HS ) 1 



(9) 



where (3 — 1/kgT and r\ is the packing fraction (77 — 
wp/6). 

All the required information to determine g(r) (a and 
pressure as functions of p) is now provided by the mini- 
mization of any of the well known and accurate density 
functionals for the Helmholtz free energy of HS solid. 
With these data, the three proposed sum rules form a 
non-linear system of equations which is solved to find A, 
ct\ and n at each p. Simultaneously, using a and Eqn.(|), 
the successive peaks of g(r) are obtained. Very recent 
studies have shown that the equation of state of the HS 
solid deduced from different functional approaches agree 
quite well with simulation results over all the density 
range ||. For the following calculation, we use the gen- 
eralized effective liquid approximation (GELA) as it 



gives the best overall behaviour. However, there is not 
significant differences if any other functional approach is 
used. If we compare the most recent Monte Carlo sim- 
ulations by Choi et al. with the predictions of the 
present theory for g(r), the agreement is excellent over 
all densities and specially impressive at high densities. 
Figures 1 and 2 show g(r) for two significant densities: 
77 = 0.52, the lowest density below melting [r\ m « 0.54) 
with available simulation data, and r\ — 0.71, near close 
packing (ij cp w 0.74), respectively. There are some dif- 
ferences between theoretical predictions and simulation 
results. Nevertheless, they are quite small and can only 
be appreciated more easily at the lowest densities. Let us 
first pay attention on g^ (r). The value at contact differs 
from that of simulation. It is a direct consequence, via 
virial equation, of the approximate theoretical pressure. 
If the exact pressure (from simulation) is used into the 
theory, the agreement with simulation would be almost 
complete (see figure 1) confirming the goodness of the 
theory. The rest of the peaks of g(r), which only depends 
on the parameter a (a = 113 for r\ = 0.54 and a = 10094 
for i] = 0.71), also agree quite well with the simulation 
results. Introducing again the exact a values (estimated 
from simulation: a = 91 for r\ = 0.54 and a = 7659 for 
r) = 0.71) the agreement is excellent. The Gaussian pa- 
rameter otx is approximately half than those of the rest 
of the peaks (ai = 50 for r\ = 0.54 and a± = 5405 for 
i] = 0.71). However, because the cutoff the width of this 
first peak is similar to the rest. The excellent agreement 
of the first peak showed in the inset of the figures would 
give an estimation of otx from simulation practically equal 
to that predicted by the theory. More interesting is the 
r\ parameter, it corresponds to the position where the 
first peak has its practically maximum value. Notice, 
however, that at very low densities the real maximun 
is located at contact (see figure 1) and r% would be the 
maximun if the first peak is analytically extrapolated be- 
low r = 1. The important point is that simulation data 
of this maximum has been reported (at very low densi- 
ties the extrapolation has been also reported). Figure 3 
shows the these data for different densities together with 
the theoretical predictions. 

According with the decreasing of the lattice parameter, 
the width of the first peak decreases monotonously with 
p {a,\ increases monotonously with p). 

r% must be always smaller that the < r > except at the 
close-packing limit where both coincide. As p decreases 
from this limit the pressure decreases rapidily from in- 
finity and therefore the value at contact. Meanwhile the 
lattice parameter hardly changes and the peak width still 
remains quite sharp, under these circumstances the only 
way to keep on with the normalization is increasing r%. 
However, at low p, the pressure does not change too much 
with p and the peak becomes duller as p decreases, r\ 
must recede to maintain the normalization. We remark 
this because the overall agreement of the peaks does not 



necessarily implies the nice agreement of the behaviour 
of r\ with the mean density. This proves the suitable 
physical description of our theoretical approach. The pa- 
rameter A is a simple factor to adjust the normalization 
or to adjust the value at contact. 

A major consequence of the theoretical knowledge of 
g(r) is the possibility to develop and use proper pertur- 
bation theories for solids where the perturbative term 
uses the appropriate correlation function of the reference 
HS system instead of a the correlation of a fluid at same 
effective density ||. Working along this direction is in 
progress. Moreover, these theories provide a way to de- 
termine g(r) for any simple system in the same way as 
in theory of simple liquids where the HS system is used 
as reference system Q . An alternative and fresh method 
is to apply directly the present approach. It can be pro- 
posed a parametric form of g(r) which includes all the rel- 
evant physics. Extending the approach to systems with 
significant compressibility, the two parameters of each 
peak, a and Ri, should differ from their homologous of 
go(r) in an amount which should decrease quadratically 
with the distance. The normalization and equilibrium 
equations should be enough to determine this decay and 
to reasonably describe g(r). Notice that the energy equa- 
tion is another sum rule which can be applied to these 
systems. The presence of defects, vacancies and inter- 
sticials, would change the normalization of the peaks, in 
both g(r) and <?o(r) which also would depend on the dis- 
tance. For this contribution, however, one should expect 
an exponential decay as the defects would behave as a 
kind of fluid inside the solid which induce a short-range 
correlation into g(r). 
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FIG. 1. g(r) at r\ = 0.52. Solid line is the prediction of the 
present theory using theoretical data from GELA functional 
approach. Dotted line corresponds to the theoretical predic- 
tions using data from a hypothetical exact theory. Triangles 
are Monte Carlo results from Choi et al. The inset shows 
details of the first peak gi(r). 



FIG. 2. As Figure 1. but for r\ = 0.71. 

FIG. 3. Parameter n (lower curve) and < r > for the first 
neighbour (upper curve) as a function of the mean density 
predicted by the theory. Triangles from simulation data. 



